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ABSTRACT 

The multiplicities of stars, and some other properties, were collected recently by Eggle- 
ton & Tokovinin, for the set of 4559 stars with Hipparcos magnitude brighter than 
6.0 (4558 excluding the Sun). In this paper I give a numerical recipe for constructing, 
by a Monte Carlo technique, a theoretical ensemble of multiple stars that resembles 
the observed sample. Only multiplicities up to 8 are allowed; the observed set con- 
tains only multiplicities up to 7. In addition, recipes are suggested for dealing with 
the selection effects and observational uncertainties that attend the determination of 
multiplicity. These recipes imply, for example, that to achieve the observed average 
multiplicity of 1.53, it would be necessary to suppose that the real population has an 
average multiplicity slightly over 2.0. 

This numerical model may be useful for (a) comparison with the results of star and 
star cluster formation theory, (b) population synthesis that does not ignore multiplicity 
above 2, and (c) initial conditions for dynamical cluster simulations. 
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1 INTRODUCTION 

Recently Eggleton & Tokovinin (2008; ET08) counted the 
stellar systems with multiplicity from 1 to 7 in the set of stars 
brighter than mag. 6 on the Hipparcos scale. There are 4558 
such sytems in total (excluding the Sun), with 2716, 1438, 
285, 86, 20, 11 and 2 having multiplicity 1 to 7. We refer 
to this as the Hipparcos-Bright Multiple System Catalogue 
(HBMSC). It is, of course, provisional: a few new compo- 
nents are found every year. The present paper proposes a 
simple model distribution that will generate the observed 
population, including the multiplicities. Selection effects are 
modeled, so that allowance can be made for the fact that 
some potential visual binaries or sub-binaries may not be re- 
solvable (depending partly on distance) , and some potential 
spectroscopic binaries may be undetectable because of unfa- 
vorable luminosity or mass ratios, or broad lines as in many 
OB stars. I estimate that to obtain the observed multiplicity 
fraction of 1.53 (ET08), i.e. an average of 1.53 components 
per system, I need an actual multiplicity fraction of about 
2.0 or more. I attempt to fit the distribution of periods (and 
sub-periods), as well as of multiplicity and period. 

Statistics on multiplicity are useful for at least three 
reasons. 

(a) They are a statement about the final product of star for- 
mation, and thus something of a check on theoretical models 
of this difficult process. 

(b) There are some stellar evolutionary scenarios that are 
importantly modified by, or even wholly dependent on, the 
existence of triple companions; thus we need to know some- 
thing of the probability of triple (or higher multiple) com- 



panions. 

(c) They represent a starting point for N-body dynam- 
ical simulations of clusters, and more generally for any 
population-synthesis project. 

Note however that the HBMSC can provide statisics only 
for systems with mass above about 1 Mq , since very few 
systems of lower mass are included among the bright stars. 
Although many low-mass multiples have been recognised, 
there does not exist a comparable wealth of data on, say, 
the 5000 nearest stars. 

In relation to point (a), it can be questioned how rep- 
resentative are the systems in the general Galactic field of 
what is produced directly in Star Formation Regions (SFRs) . 
This is not the place to attempt a detailed discussion of star 
formation, but there seems no doubt that the great major- 
ity of stars form in localised, relatively dense, SFRs, and are 
then released or scattered into the Galactic field. Since only 
about 1% of stars are currently in clusters, it follows that the 
field stars are the great bulk of what is produced in SFRs, 
and so it is not unreasonable to see them as representative 
of those products. Some of the systems have probably been 
modified, by N-body encounters, during both the formation 
process and the later escape process. One might attempt to 
distinguish between the formation of stars and the forma- 
tion of systems, but although a semantic distinction might 
be made I doubt if any clear physical distinction can be 
made. The widest systems that can be found in the field are 
usually much too wide to have had a chance to form in the 
relatively compact SFRs where they were born, but presum- 
ably formed during the process of evaporation that put them 
into the field. Those field systems that are not very wide can 
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reasonably be supposed to be much the same as when they 
formed in the SFR. Numerical attempts to model star for- 
mation in SFRs will presumably have to address the N-body 
interactions and the dissipation into the field as well as the 
gas-condensation stage, and so it is not unreasonable that 
they should be constrained to some extent by the statistics 
of multiplicity in field stars. 

These remarks also relate to point (c). The distribution 
of multiples in a cluster will certainly differ from that of 
the field, by missing the widest field multiples. But if one 
starts an N-body cluster with the multiplicity distribution 
suggested below, the widest systems will be very quickly 
broken up. It is not unreasonable to suppose that the closer 
systems, which do not get broken up quickly, might be rep- 
resentative of what should be expected in a young cluster. 

In Section 2 I describe the Monte Carlo procedure I use, 
applied to a model in which stellar evolution, but only so far 
of single stars, is taken from a table of evolved stellar mod- 
els. In Section 3 the procedure is generalised from single to 
multiple stars; a procedure called Create constructs a the- 
oretical population of multiple systems. Although the issue 
of interactive binary and triple systems is briefly discussed, 
it is left for future work. In Section 4 I describe a 'theoreti- 
cal observatory', a procedure called Observe, which looks at 
the multiple systems generated by Create and decides which 
systems and subsystems should be recognisable as visual or 
spectroscopic or eclipsing binaries, and which not. Thus a 
'raw' model from Create is turned into a model that can 
in principle be compared with what is actually seen, in the 
HBMSC. In Section 5 a comparison is made, and I discuss 
some applications of the model to the problems mentioned 
in the second paragraph. 

I would like to emphasise that I am not claiming to 
find a definitive model of multiplicity, but only a reasonable 
model. For instance, a definitive model would need many 
more mass ratios than are currently well-known, and in par- 
ticular rather extreme mass ratios, such as are very difficult 
to determine reliably at present. In addition, there is an ob- 
vious degeneracy between the Create stage and the Observe 
stage, since one can imagine that many more binaries (and 
higher multiples) are created that have parameters that are 
hard to observe, and so are excluded at the Observe stage. 

2 METHOD 

We would like to model the statistics with some set of (cumu- 
lative) distribution functions. Suppose for the sake of argu- 
ment that all stars are single stars. Each star is determined 
by two parameters, mass (m) and age (t). I ignore metalic- 
ity for the time being - and in fact throughout this paper - 
although it would not be difficult in principle to include it. 
For determining which stars are visible above some apparent 
magnitude, a third parameter, distance (d), has also to be 
assigned. The cumulative distributions of these parameters 
will be some functions of the three variables. There is no a 
priori reason to suppose separability, i.e. that the distribu- 
tion is the product of three elementary distributions, each in 
one variable only. We select three random numbers X, Y and 
Z, uniform in the interval (0, f ), and then generate m from 
some global (cumulative) IMF, m = m(X), t from an age 
distribution t = t(Y, m), and finally d from a distance distri- 
bution d — d(Z,m,t). The order is immaterial, although of 



course the actual functional forms will depend on the order 
selected. The distribution over distance will have to fall off 
at large distances fast enough to avoid Olbers' paradox; the 
formulation given below - equation (10) or its inverse equa- 
tions (12) and (13) - is roughly spherical within ~ 100 pc, 
roughly disk-like to ~ 1000 pc and becomes negligible be- 
yond ~ 3kpc. 

The number of stars brighter than a star of luminosity 
Lo at distance do, i.e. that have 



(1) 



d? ~ d% ' 

out of a total population of N* stars is an integral 

N = N* [ H{L(X,Y) - Lo^jjP-} dX dY dZ , (2) 

J unit cube 

where H is the Heaviside step function. The integral can be 
rewritten as 

N = N* Z v (X,Y)dXdY , (3) 

J unit square 

where Z v is the fraction of stars at distance less than d v , 
and d v is the distance out to which a particular star as a 
funtion of X, Y is visible. This distance d v is given by 



J L(X,Y) 

«0 7 



(4) 



Given that stars are discrete rather than continuous, we 
replace the integral by a Monte Carlo approximation. Sup- 
pose we select n 2 points randomly in the X, Y unit square. 
Then we approximate the integral (3) by 



i=i 



(5) 



It is clear that only stars that come from a rather small 
fraction of the X, Y plane will populate the night sky visibly 
to the naked eye: we know from experience that only stars 
more massive than about 0.8 Mq, i.e. ;€ 10% of stars from a 
reasonable IMF, can be bright enough, and the most massive 
and brightest stars must be very young, e.g. age less than 
~ 10 7 yrs at masses ^ 10 Mq . We can expect that a rather 
small fraction of the n systems will contribute at all, but 
some of these will contribute rather heavily. We therefore 
use 'importance sampling', i.e. a transformation that maps 
the unit square in X, Y from the unit square in say X' , Y' so 
that the region near X = 1, Y = 0, where most visible stars 
come from, is sampled more densely. Following Eggleton et 
al. (1989; EFT89) I use 



X = 1 - (1 - X'j 2 



l-(i-y') 



'n(I-X') 



(6) 



In equation (5) we have to insert an extra factor, the Jaco- 
bian of the transformation: 



d(X',Y') 



(7) 



Equation (7) gives a fractional number of stars at each 
sampling point. We quantise this in the following way. At the 
ith element the value Nt,Z v J/n 2 in equation (7), where J is 
the Jacobian, is some fractional number, say k + f, where k 
is an integer and < / < 1. Then select a further random 
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number, uniform in [0, 1). If this is less than / then k + 1 
'cloned' stars are allowed, but if greater, then k clones. A 
value for n of 200 was found to give a reasonable compromise 
between speed and detail. 
The formula 

/ X \ 0,85 

m = 0.3 [y^y) , (8) 

gives a distribution which, in terms of log m, is peaked at 
X = 0.5, and falls off with a moderately Salpeter-like slope 
of 2.23 for m^O.8M0. The median mass, at X = 0.5, is 
m = 0.3 Mq. For age we use the simplest distribution 

t = Y, (9) 

t being the age in units of 10 10 yrs. Although for a total 
Galactic population such a uniform rate of star production 
is unlikely to be realistic, most of the bright stars are less 
than a tenth of the Galaxy's age, and a uniform birth rate 
is not so unreasonable in this span. 

Formulae (8) and (9) can be seen as 'plucked out of thin 
air'. However a distribution like (8) has the desirable qual- 
ity of resembling the Salpeter distribution at large masses, 
but turning over in a reasonable way at low masses, be- 
low 0.3 Mq. The value of the turnover mass and of the slope 
below that mass are not determinable from the HBMSC, be- 
cause the latter includes hardly any stars below ~ IMq. 
The slope at the high-mass end is in principle determinable, 
say by doing some least-squares fitting; but my initial value 
of 0.75 seemed to produce too many high-mass relative to 
intermediate-mass stars, and a second guess at 0.85 seemed 
to be about right, as quantified by a \ 2 test given near the 
end of this Section. The assumption in (9) of a reasonably 
uniform star formation rate over the last Gyr is similarly 
bound to be incorrect in detail (i.e. on a time scale of Myrs, 
and a lengthscale of tens of parsecs), and yet is a fairly rea- 
sonable model on longer timescales and lengthscales. 

The cumulative distribution Z(d) is approximated by 



fto ft? 



1 

Z 

so that 



(10) 



Z txd 3 if d £ ft < hi 

oc d 2 if ho&d&hi 

~ 1 if dS; hi . 
We need the inverse of this, d(Z). If 



3ho f 1 3(1 - Z) 
2hi ' 



(11) 



(12) 



then 



3ft, 

x 



,1 -1 \ -c 

— cos(- cos X) it 



< 1 



cosh(— cosh 1 x) 
x 3 



if 



> 1 



(13) 



A refinement that we make to this d(Z) model is that 
we take fto, which is essentially the scale height of the disc, 
to be dependent on age: 

h (t) = ht 3 , h = 200 , (14) 



with ft, fto in parsecs and t in units of 10 10 yrs, as before. 
The constant fti is taken as 1 kpc. Another refinement is 
that reddening, at 1 mag. per kpc, is included. 

Once again formulae and coefficients are apparently 
plucked out of thin air, but are based nevertheless on a 
general picture of the Solar neighbourhood that I believe 
is widely accepted: that the distribution is roughly spherical 
for old stars but more disc-like for young stars. I have sim- 
ply adopted what I believe is about the simplest formulation 
that describes this. I make no attempt to improve the fit by 
varying some coefficients. 



a) Theoretical HRD 



O 
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Fig 1 - Hertzsprung-Russell diagram for 171 stellar 
models with initial masses in the range 0.1 — 63 Mq. Evo- 
lution was terminated when it became rapid, approaching 
a final compact remnant; except that an upper age limit of 
10 Gyr was also imposed. 

I used my stellar-evolution code (Eggieton 1971, 1972, 
Eggleton et al. 1973, Eggieton 2006; Pols et al. 1995) to 
evolve 171 single-star models with (log) masses of 



log m = - 1.0 (0.02) 0.0 (0.01) 0.6 (0.02) 1. 



(15) 



All had a roughly solar composition: X — 0.70, Y = 
0.28, Z = 0.02. They were evolved from the Zero-Age Main 
Sequence (ZAMS) either until 10 10 yrs old or until fairly near 
the end of their lives. For the lowest masses this meant 
very little evolution in practice. For intermediate masses 
I stopped the evolution either at the very rapid phase of 
evolution towards the white-dwarf region, or shortly after 
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the onset of carbon burning in the core. A rate of mass 
loss by stellar wind was adopted (Eggleton 2006), so that 
intermediate-mass stars less massive originally than 4.2 Mq 
would become white dwarfs less massive than 1.05 Mq. All 
estimates of mass-loss rates for red supergiants are very ten- 
tative; the one used here gave relatively little mass loss in 
the (initial) mass range 6—20 Mq, but proportionately much 
increased mass loss at still higher masses. 

For every tenth timestep of each evolutionary sequence, 
the following six numbers were stored: t, m(t), log L, log R, 
log i?max , and I cv . -R max is the maximal radius in the evo- 
lution prior to age t, useful for knowing whether binary in- 
teraction with a fairly close companion will have occurred, 
and lev is an integer that describes crudely the evolution- 
ary state, from 1 for an MS star to 10 for a neutron star or 
black hole. Fig. 1 shows a theoretical Hertzsprung-Russell 
diagram (HRD) for this collection of stars. To avoid some 
interpolation which would often be of a highly non-linear 
nature, when the Monte Carlo procedure of equation (8) se- 
lects an initial mass I replace it by the nearest value in the 
sequence (15). 

Table 1 presents some results. The distributions (8) for 
mass, (9) for age and (10) for distance were used, and only 
N*, the total population, was varied to get about the right 
number (to much better than 1%) of visible stars. The distri- 
bution over spectral type (Table 1, row labeled 'single'), to 
some extent equivalent to the distribution over mass at least 
for types F and earlier, is not very satisfactory compared to 
the observed distribution. The ratio computed/observed is 
nearly 3 for O stars, and ranges between 1.4 and 0.7 for the 
remaining spectral types. 

Although we might try to ameliorate this by playing 
with the IMF, it is not difficult to see that much of the 
disagreement is really a consequence of the fact that the 
distributions of mass, age and distance are not independent 
of each other, as is implicit in equations (8) - (10); although 
equation (10) for the distance dependence does include an 
age dependent term via equation (14). The regions that form 
massive stars are very clumpy, with several hundred parsecs 
between them, and the Sun appears to be located within a 
'bubble' (Cox & Reynolds 1987) of 100 - 250 pc in which 
no formation of massive stars has taken place recently, i.e. 
within about 75 Myr, the lifetime of a 6 Mq star. O stars 
have formed much more recently than that within say the 
Orion SFR (Star Forming Region) at about 450 pc, and are 
easily above H p — 6, but there is no O star nearer than 
about 250 pc. Thus one way to improve the agreement at O 
and early B is to exclude stars above 6 Mq that are nearer 
than 250 pc, and that was used in the line of Table 5 labeled 
'OB down'. This is equivalent to saying that the mass, age 
and distance distributions are not in fact disjoint; but it is 
numerically easiest simply to omit systems that have mass 
above a threshold and distance below another threshold. 

We might try some more alterations, particularly of the 
IMF, to improve agreement further, but in practice much 
of the remaining disagreement disappears when we include 
multiple systems, in the next Section. The row labeled 'mul- 
tiple' in Table 1 is essentially the same model as the 'OB 
down' model, but with multiplicity included as will be dis- 
cussed in the next Section. The agreement is either slightly 
or in some cases considerably better. The overall improve- 
ment is at least partly because massive systems, as we will 



see shortly, appear to be more highly multiple than less mas- 
sive systems, and so some O and eB systems are broken 
down to later types. A \ 2 test comparing the four distribu- 
tions over spectral type listed in Table 1 gives, of course, 
very poor agreement between the observed distribution and 
either of the first two theoretical distributions, 'single' and 
'OB down'. But we get a \ 2 of 11.2 when comparing the 
'multiple' with the 'observed' distribution, which is accept- 
able at the 10% level. No doubt this could be improved by 
tweaking some coefficients, but this does not seem necessary. 

For each of the 'single', 'OB down' and 'multiple' lines 
in Table 1, ten Monte Carlo models were run differing only in 
the random numbers generated. What is listed is the average 
of these ten, and the rms variation. A few preliminary runs 
were done in order to estimate the value of the normalisation 
constant N* of equation (7) that would give a final count 
(averaged over the ten runs) of approximately the number 
4558 observed. 



3 A MODEL OF MULTIPLICITY 

When deciding on the multiplicity of a computed stellar 
system, I continue to use random numbers, thinking in terms 
of a process of successive bifurcation; although for practical 
reasons I allow no more than 3 successive bifurcations so 
that only multiplicity up to octuple is considered. Note how- 
ever that this is not a claim that real multiples are formed 
by successive fragmetation during a pre-main-sequence con- 
traction phase. It is simply an acknowledgment of the well- 
established fact that the great majority of multiple systems 
are highly hierarchical (Evans 1968), and such systems by 
definition consist of binaries nested inside wider binaries. 
A rather small proportion of systems is non-hierarchical, in 
the sense that 3 or more components are found at roughly 
equal distances from each other. Seventeen such systems 
were identified in the HBMSC of ET08, but these are among 
the least certain. Their inclusion or exclusion will make little 
difference to the following. 

We have already selected a mass m, which we now in- 
terpret as the total mass, and the age t. When considering 
whether to bifurcate or not, we choose two more random 
numbers uniformly probable in (0,1), U and V say, to deter- 
mine (i) the mass ratio Q (with a finite probability of zero, 
i.e. no bifurcation) and (ii) orbital period P. Of course, both 
P and Q are irrelevant if the choice of U is found to imply 
no bifurcation after all. P and Q might, indeed probably 
will, depend on m as well as U, V , and also on the hierar- 
chical depth I, i.e. on whether we are creating a binary, a 
sub-binary or a sub-sub-binary. We allow only I = 0, 1 or 2. 
At each potential bifurcation we choose a fresh U and V, so 
that we may need as many as 7 U, V pairs, although with 
some probability that there may be no bifurcation at all, 
and only one U is used. 

First we determine whether the current hierarchical 
level I bifurcates or not. I specify (Table 2) a tabular func- 
tion Uo(m,l), I being the hierarchical depth. Then U > Uo 
implies no bifurcation, and U < Uo implies bifurcation with 
a mass ratio that we can suppose depends in some manner 
on (Uo — U)/Uo, a ratio which is always between zero and 
unity, and with uniform probability. Actually the entries in 
Table 2 are treated as functions of a variable j related to 
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Table 1. Distribution over spectral type. 



O eB IB A F G KM tot. 

observed 38 401 653 928 578 587 1042 331 4558 

single 97 620 897 886 468 492 881 225 4565 

rms 16 44 46 65 25 28 79 87 138 

OB down 83 440 929 915 491 509 923 260 4550 

rms 18 42 60 41 42 24 95 102 123 

multiple 47 402 699 917 581 529 1041 335 4550 

rms 9 45 55 39 34 35 101 56 148 

eB stands for early B (B0 - B3.5), and IB for later B. 
'Single' refers to the single-star model of Section 2. 

'OB down' is a model, described at the end of Section 2, where stars over 6Mg and nearer than 
250 pc are discounted. 

'Multiple' is the same as 'OB down', but with multiplicities up to 8 allowed as in Section 3. 



Table 2. Bifurcation Probability Uo as a Function of Mass and Hi- 
erarchical Level. 



m = 





.01 


.09 


.32 
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3.2 


11 


32 


oo 


3 = 


1 


2 


3 


4 


5 


6 


7 


8 


9 


I = 


0.40 


0.40 


0.40 


0.40 


0.50 


0.75 


0.88 


0.94 


0.96 


1 


0.18 


0.18 


0.18 


0.18 


0.18 


0.18 


0.20 


0.60 


0.80 


2 


0.00 


0.00 


0.00 


0.00 


0.00 


0.20 


0.33 


0.82 


0.90 


3 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 



The values for m^l are notional, because very few systems of this low a 
total mass appear in the HBMSC. However they are reasonably consistent 
with what is known about the nearest 66 systems. 



total mass m by 



instance, 



5 + 2 log 



100m N 



100 + m J 



(16) 



The variable j ranges from 1 to 9 as m ranges from zero 
to infinity, and I interpolate linearly in the Table for non- 
integer values of j. 

Then for I = 0, and assuming that the random U does 
imply bifurcation, a new random number V is selected in 
(0,1). I take the period to be 



10° 



V 2 



(i-vy 



(17) 



in days. This expression gives a distribution of log P peaking 
at about 10 J d, dropping rather slowly towards large P and 
also dropping, a little more rapidly, towards small P. The 
lowest ten percentile is at ~ 1.3.10 3 d, and the highest at 
2.5.10 7 d. For I > 0, I adopt 



P = 0.2Pi_ilO" 



(18) 



The median period ratio is 10 -3 ' 2 . These determinations 
of outer or inner periods may well be expected to depend 
also on all previous choices, in particular on m and I. They 
might depend on t as well, if we are talking about a cluster 
in which orbits can soften or harden as time advances; but 
for the time being I do not allow such complexity. 

The mass ratio Q can be determined from the previ- 
ous U and corresponding bifurcation probability Uo by, for 



Q = 



( Uo-U 
\ Uo 



a = f(m,l,Pi) 



(19) 



whether I is 0, 1 or 2. Values of Q less than 0.01 are replaced 
by 0.01. Because, as noted above, (Uo — U)/Uo is uniformly 
distributed in (0,1), distribution (19) is just a power law. Ob- 
viously with enough information it will be better to have a 
more complicated distribution, but several attemped deter- 
minations of the mass-ratio distribution have been presented 
as power laws. For example, Kouwenhoven (2006) obtains, 
for near- infrared companions to IB/ A stars in the Sco OB2 
association, a differential distribution oc Q . This cor- 
responds to a ~ 1.5 in equation (19). Tokovinin (2000), 
following Lucy & Ricco (1979), finds a marked preference 
for 'twins', i.e near-equal-mass pairs, at fairly short periods 
( ^ 25 d). This would require a low power like a ~ 0.15 if one 
wanted 50% of systems to have Q > 0.9. For the present I 
assume, following EFT89, that a = 0.8 for all periods above 
25 d, and for 75% of those with shorter periods, but for the 
remaining 25% of the latter I take 



Q = 0.9 + 0.09 



( Uo-U 
V Uo 



(20) 



One would like to do the same analysis for the 4500 
nearest stars, but although information on these stars, par- 
ticularly on L and T dwarfs, has increased by leaps and 
bounds over the last decade, there is nothing like as com- 
plete a study of their multiplicity as already exists for the 
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~ 4500 brightest stars. However, for the 66 nearest stars, 
including the Sun, multiplicities 1, 2, 3 and 4 have frequency 
42, 16, 8 and 0. These are roughly the values we expect using 
the values for Uo(m, I) in the columns of Table 2 appropriate 
to m 1. 

Obviously any bifurcation will have to be reversed, i.e. 
ignored, if the period selected turns out to be too small 
for the components to exist separately. This may be seen to 
beg the question of contact binaries. I have argued elsewhere 
(Eggleton 1976, 1996) that contact binaries are not 'born' in 
contact, but their periods evolve towards contact as a result 
of one or both of the following processes: magnetic braking 
and Kozai cycles, both assisted by tidal friction. I return to 
this point in Section 5. The question is slightly academic 
in the present context, as there are only three contact bi- 
naries in the HBMSC; but mergers may mean that there 
are some former contact binaries as well. However it is ev- 
idently necessary to place some lower limit on P based on 
the sum of the radii of two ZAMS stars relative to their sep- 
aration. It is also necessary to have some upper limit. I take 
the upper limit to be 10 10 d (27Myr) for hierarchical depth 
I = 0, and, according to equation (19), 0.2P(_i for hierarchi- 
cal depth I > 0. The former comes loosely from the consid- 
eration that wider orbits would be disrupted by independent 
neighbouring stars, and the latter from the dynamical sta- 
bility of coplanar triples (Harrington 1975). All orbits are 
assumed circular for the present. It would be easy to ran- 
domise the eccentricities if their cumulative distribution is 
postulated; but this would introduce the issue of modeling 
tidal friction to vary the eccentricity with time, and for the 
present I want to adhere to simplicity. The lower limit for Pi 
comes, as mentioned before, from the condition that a star 
should not fill its Roche lobe on the ZAMS. 

A set of systems formed by this procedure will be 
called a 'Raw Theoretical Bright Multiple Star Catalogue' 
or RTBMSC; in the next Section I describe the 'Theoreti- 
cal Observatory' which decides which systems, subsystems 
and subsubsystems should be recognisable, and which not. 
The result of that process will be a Theoretically Observed 
Bright Multiple Star Catalogue, or TOBMSC. Some sys- 
tems of particularly high luminosity appear several times 
over in the RTBMSC, according to the 'cloning' recipe be- 
tween equations (7) and (8). Create assigns a random incli- 
nation to each pair, at each hierarchical level, and does this 
independently for each system, including each clone. Such 
random inclinations should in fact influence the dynamical 
stability, but I do not include such an effect yet. For the 
present, they influence only the likelihood that a subsystem 
will be detectable as either a spectroscopic or eclipsing bi- 
nary by the Theoretical Observatory (Section 4). 

Equations (17) to (20) are simpler than might have been 
expected. We could expect that the probability distribution 
of period might depend on the mass, and that similarly the 
distribution of mass ratio might depend on mass and the pe- 
riod. My prescription (19), (20) for Q is period- dependent, 
as described in the text, but the effect of this period depen- 
dence is fairly modest in view of the modest numbers that 
have P < 25 d. 

The computer code Create contained in principle 45 
adjustable constants, although only 12 of the 27 constants 
in Table 2 (corresponding to the larger elements in the rows 
/ = 0, 1 and 2) were in practice adjusted. In addition, only 



6 of the 18 remaining constants were varied, from guessed 
initial values, to try and improve the agreement. The 12 not 
varied include, for instance, almost all of the coefficients and 
exponents in equations (8) - (14) and (17) - (20). In most 
cases it was possible to estimate from an eyeball analysis 
of the HBMSC what values would be plausible, and these 
values did indeed appear to work without further tweaking. 
It would be only a slight exaggeration to say that the coeffi- 
cient and two exponents of equation (17), for example, were 
'plucked out of thin air'. But they were based on an eyeball 
analysis of the HBMSC section of Table 6, along with the 
suspicion that the rather marked shortage of systems there 
with periods of 1 - 100 yr at early spectral types might be 
an effect of observational selection, as is indeed found in the 
next Section. 

Once an n-tuple system has been created, the various 
components can be evolved, all according to the same pre- 
scription as described in Section 2. I do not (yet) put in a 
binary interaction, such as Roche-lobe overflow (RLOF). It 
is feasible to put in back-of-the-envelope prescriptions (Hur- 
ley et al. 2002), but I prefer, in a future version, to utilise 
the capacity of my stellar evolution code to follow RLOF in 
some circumstances (Nelson & Eggleton 2001, de Mink, Pols 
& Hilditch 2007) . But a substantial proportion of potentially 
interacting systems are quite difficult to follow numerically, 
because the RLOF can become very rapid (Paczyriski & 
Sienkiewicz 1972). In addition winds, including those linked 
magnetically to a component and capable of draining angu- 
lar momentum from the component, and via tidal friction 
from the binary, can affect the orbit. Such winds are in- 
cluded in my binary evolution code, which can work in a 
mode where both components, and the orbital period, and 
in addition the eccentricity, which is also modified by tidal 
friction, can be solved for in a single implicit time step. This 
process even extends to models in which both stars over- 
fill their Roche lobes, and influence each other's structure 
via internal luminosity transfer as is frequently hypothesised 
for contact binaries (Lucy 1968, 1976, Yakut & Eggleton 
2005); but this mechanism still contains some very uncer- 
tain physics, and it is also unusually expensive to compute 
since evolution progesses in cycles on a thermal timescale, 
while lasting for something like a nuclear timescale. Thus 
the number of timesteps is larger than usual by a factor of 
several hundred. 

The worst cases of rapid mass transfer, usually hy- 
pothesised to lead to common-envelope evolution (Paczyriski 
1976), and to either very close white dwarf + MS star bina- 
ries or else to a merger, are not dealt with yet, even by the 
most sophisticated version of my evolution code. In practice, 
these are likely to be the majority of RLOF cases, and so 
for the present I simply annotate those systems in which one 
component is, or has been at some time in the past, larger 
than its Roche lobe. 

Quite a few massive systems will have produced super- 
novae by now, even though one or two OB stars (and even 
AF stars) remain in the system. This is partly why we need 
multiplicity significantly higher in O stars than in much 
later stars. Although we expect that such neutron stars 
should have ejected themselves from their natal systems by 
means of asymmetric mass ejection (Shklovskii 1970, Lyne 
& Lorimer 1994, Hansen & Phinney 1997), I choose in the 
present work to keep them in the system so that, for ex- 
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Table 3. Velocity and Temperature Thresholds as functions of Spectral Type. 

Type O eB IB A F G KM 

V cr (km/s) 30.0 25.0 20.0 12.0 6.0 2.0 2.0 3.0 

T eS (K) 31000 15500 10200 7400 5900 4880 4000 

The first row of numbers is the minimal radial velocity amplitude considered measurable at that spectral 
type. The second row is the minimal effective temperature considered belonging to the spectral type. Note 
that the G/K temperature boundary (4880K) has to be appropriate to giants rather than to dwarfs. 



Table 4. Multiplicity Frequency. 



Sample 


n = 1 


2 


3 


4 


5 


6 


7 


8 


tot. 


av. 


HBMSC 


2716 


1438 


285 


86 


20 


11 


2 





4558 


1.53 


RTBMSC 


1459 


2179 


517 


202 


101 


44 


30 


18 


4550 


2.04 


rms 


63 


111 


37 


28 


10 


11 


13 


7 


148 


0.03 


TOBMSC 


2649 


1472 


302 


85 


29 


9 


3 





4550 


1.55 


rms 


96 


85 


37 


17 


5 


3 


3 





148 


0.02 



The last column is the average multiplicity or 'companion star fraction (CSF)', defined as 

J2 n nN n/ J2n Nn - Ten models 

were generated by Create, averaged to give the line RTBMSC, 
then theoretically observed by Observe, and finally averaged to give the line TOBMSC. 



ample, I can estimate from the RTBMSC how many neu- 
tron stars should have been produced by the systems of the 
HBMSC: 220 ± 24. Many more systems will have produced 
white dwarfs, and these are rather more likely to have re- 
mained in their original systems, although some, perhaps 
the majority, may have merged with a companion star. The 
number of white dwarf companions, seen or not, and merged 
or not, is 371±60 in the RTBMSC. Only 20 are known in the 
HBMSC, which suggests that many WD companions have 
not so far been recognised as such. 

4 THE THEORETICAL OBSERVATORY 

I use the code described in the previous Section and called 
Create to create a 'raw' catalogue of theoretical bright 
multiple systems, the RTBMSC. Then I use a code called 
Observe which reads in the output catalogue from Create 
and decides which components would actually be identi- 
fiable, either as visual systems, spectroscopic systems or 
eclipsing systems. This works on very simple principles, out- 
lined below, and converts the kind of statistics seen in the 
middle sections of Tables 4, 5 and 6 into those of the bottom 
sections. Let us call this the 'Theoretically Observed BMSC 
or TOBMSC. 

For visual binaries to be resolvable we require that one 
or other of 

AH P < 2.5(p-0.l") , (21) 

AH P < 4.5 log(p/0.05") , (22) 

be satisfied, where p is the separation in arcsecs and AH P 
is the modulus of the difference in Hipparcos magnitudes. 
A further condition is that the fainter (sub)component 
should be brighter than 14th magnitude. Criterion (21) is 
roughly appropriate to normal visual measurements, and 
the much stronger criterion (22) is roughly appropriate for 
high-resolution measurements with speckle or adaptive op- 
tics (Sterzik & Tokovinin 2002). For the present I assume 



simply that 50% of stars have been examined by one pro- 
cess and 50% by the other, chosen at random. 

For a spectroscopic binary (SB) to be identifiable I as- 
sume that the radial velocity amplitude of the brighter com- 
ponent, Ki, is above a threshold value V CI that is given in 
Table 3, and further that the period is less than ~ 50 yr: 

> V„ , P < 18000 d . (23) 

The threshold V cr depends on spectral type. Theoretical SBs 
can be single or double lined, depending on the luminosity 
ratio. The value of K\ includes an inclination which was 
assigned randomly to each sub-binary in the Create stage. 
This implies the assumption that inclinations at different 
hierarchical depths are uncorrelated. It is hard to test this: 
Muterspaugh et al. (2006) showed that only six triples have 
unambiguous determinations of the inclination of one orbit 
to the other (as distinct from the inclination of each orbit to 
the line of sight). It is unlikely that 6 data points would give 
a clear determination of the statistics. In fact the inclinations 
were not very consistent with a complete lack of correlation, 
nor with complete correlation. Sterzik & Tokovinin (2002) 
considered 135 visual triples, and estimated the mean angle 
between orbits as 67° ± 9° . This value would be 90° if the 
two angular momenta were uncorrelated; thus there appears 
to be some correlation but not a great deal. However, for 
most of these systems there was a 180° ambiguity in the 
position angle of one or other line of nodes, which creates 
an ambiguity in the relative angle of the angular momenta; 
this ambiguity was missing in the much smaller sample of 
Muterspaugh et al. (2006). 

The inclination assigned in Create also says, along 
with the period and radii, whether a system or subsystem 
will eclipse, and further gives an estimate of the depth of 
eclipse. Provided the eclipse depth is sufficient, such eclipsers 
are added to the tally of recognised theoretical binaries. 
However all theoretical eclipsers turn out to be theoreti- 
cal SBs, whether single or double-lined. I do not at present 
count ellipsoidal variables, although they are included in the 
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Table 5. Multiplicity Frequency by Spectral Type. 
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tot. 
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HBMSC 






















O 


13 


10 


7 


4 


2 


2 








38 


2.42 


eB 


213 


125 


46 


8 


6 


1 


2 





401 


1.70 


IB 


353 


215 


61 


21 


1 


2 








653 


1.63 


A 


501 


336 


63 


17 


6 


5 








928 


1.61 


F 


304 


219 


36 


15 


4 











578 


1.61 


G 


323 


217 


33 


12 


1 


1 








587 


1.56 


K 


739 


263 


32 


8 














1042 


1.34 


M 


270 


53 


7 


1 














331 


1.21 


tot. 


2716 


1438 


285 


86 


20 


11 


2 





4558 


1.53 


RTBMSC 






















O 


7 


15 


6 


7 


6 


2 


2 


2 


47 


3.31 


eB 


76 


162 


56 


39 


38 


15 


10 


6 


402 


2.80 


IB 


247 


339 


71 


25 


9 


4 


3 


2 


699 


1.92 


A 


320 


453 


102 


27 


9 


3 


3 


1 


917 


1.88 


F 


211 


264 


67 


26 


6 


4 


2 


1 


581 


1.93 


G 


172 


269 


58 


19 


7 


3 


1 


1 


529 


1.94 


K 


351 


490 


122 


45 


17 


7 


6 


4 


1041 


1.99 


M 


77 


188 


36 


15 


9 


6 


4 


1 


335 


2.19 


tot. 


1459 


2179 


517 


202 


101 


44 


30 


18 


4550 


2.04 


TOBMSC 
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14 


15 
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6 


2 


1 








47 


2.41 


eB 


192 


126 


44 


24 


12 


3 


1 





402 


1.88 


IB 


458 


205 


28 


6 


2 











699 


1.41 


A 


561 


300 


47 


9 


1 











917 


1.46 


F 


333 


193 


41 


11 


3 


1 








581 


1.56 


G 


278 


205 


36 


8 


1 











529 


1.58 


K 


616 


325 


74 


16 


6 


2 


1 





1041 


1.54 


M 


198 


102 


23 


6 


3 


2 


1 





335 


1.56 


tot. 


2649 


1472 


302 


85 


29 


9 


3 





4550 


1.55 



Early B stars (B0 - B3.5) are called 'eB'; later B stars are called 'IB'. Wolf-Rayet stars are 
included under O; S and C stars under M. The second last column is the total number; 
the last column is the average multiplicity. 

HBMSC: the observed Hipparcos-Bright Multiple Star Catalogue (ET08). 

RTBMSC: a Raw Theoretical Bright Multiple Star Catalogue (Section 3). 

TOBMSC: a Theoretically Observed Bright Multiple Star Catalogue (Section 4). 

Each row of RTBMSC and TOBMSC is an average of 10 simulations, rounded to the 

nearest integer. Consequently the row or column for 'total' is not necessarily the exact 

total. 



HBMSC. 

5 DISCUSSION 

Tables 4, 5 and 6 present some details of the multiplicity 
model that appears to represent reasonably well the ob- 
served data of the HBMSC. Table 4 shows the distribution 
over multiplicity. The TOBMSC has somewhat too many 
triples and quintuples, somewhat too few binaries and sex- 
tuples, and roughly the right average. The last, of course, is 
not by accident; the entries in Table 2 were varied from ones 
initially guessed until the average came out about right. I 
should emphasise again that only those values in Table 2 
that relate to m (total mass) above unity are tested by the 
HBMSC. 

The top third of Tables 5 and 6 gives the HBMSC data 
from ET08. The middle third gives the equivalent data that 
emerges from the Create stage described in Section 3. The 
bottom third gives the result of processing the middle third 
through the selection effects of the Observe stage of Sec- 
tion 4. I ran ten cases differing only in the random number 



selection. Tables 4-6 present the average of these ten, and 
in some cases the rms scatter. The 10 cases gave an average 
total of 4550 ± 148 systems. Because the middle and lower 
thirds are both averages of 10 cases, the rows and columns 
labelled 'total' in Tables 5 and 6 are not necessarily the ex- 
act totals but only the approximate totals of the integers 
in the bodies of these Tables. Note that the RTBMSC in 
Table 4 produced 18 ± 7 octuple systems, but none of these 
fully survived the scrutiny of the Theoretical Observatory. 

When in Table 5 we compare the distribution over spec- 
tral type (second-last column) of the TOBMSC with the 
HBMSC we get a x 2 of 11.2, as already quoted in Section 2. 
When we compare the distribution over multiplicity (last 
row) of the TOBMSC with the HBMSC we get a x 2 of 
7.3. The first value is on the margin of significant or non- 
significant discrepancy at the 10% level, while the second 
favours no significant discrepancy. However a x 2 test is not 
really going to validate (or contradict) the model, since we 
could have the same x 2 using either (a) a model of multi- 
plicity that gives more multiples, combined with a model 
of selection effects that makes them harder to recognise, or 
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Table 6. Period Distribution in Systems and Subsystems. 



logP(yr) -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 tot. 
HBMSC 
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eB 





25 


41 


21 


14 


21 


26 


29 


33 


21 


3 





1 
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25 


53 


24 


20 


43 


62 


63 


48 


16 


8 
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A 





27 


62 


25 


46 


78 


78 


66 


47 


24 


9 








462 


F 





14 


33 


24 


39 


46 


61 


44 


26 


14 


3 


1 





305 


G 


1 


7 


9 


20 


40 


38 


49 


38 


32 


23 


6 


1 





264 


K 





4 


4 


11 


56 


35 


35 


42 


35 


26 


4 








252 


M 


1 








4 


9 


10 


4 


9 


9 


3 


3 








52 


tot. 


2 
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213 
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224 
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242 
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36 


2 


1 


1979 
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21 
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11 


15 


12 
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5 


2 


1 


1 
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eB 


1 


84 


96 


87 


89 
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99 


77 


48 


25 


10 


6 


1 
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IB 


5 


29 


42 


53 


79 
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129 
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57 


27 


9 


4 
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A 


8 


34 


52 


62 


97 
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77 


31 


12 


9 





810 


F 


8 


29 


38 


51 


64 


86 


104 


74 


50 


21 


8 


7 





539 


G 


5 


19 


33 


48 


58 


92 


87 


78 


47 


19 


7 


5 





497 


K 


4 


50 


77 


86 


122 


166 
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164 


94 


56 


25 


13 


1 


1032 


M 


1 


18 


35 


35 


36 


56 


48 


72 


60 


25 


6 


7 





399 


tot. 


32 
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206 


77 


52 


2 
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79 


85 
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10 
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42 
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44 
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3 
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178 


tot. 


2 


164 


249 


232 


233 


270 


297 


394 


275 


120 


48 


31 


1 


2315 



HBMSC: the observed Hipparcos-Bright Multiple Star Catalogue (ET08). 

RTBMSC: a Raw Theoretical Bright Multiple Star Catalogue (Section 3). 

TOBMSC: a Theoretically Observed Bright Multiple Star Catalogue (Section 4). 

Each row of RTBMSC and TOBMSC is an average of 10 simulations, rounded to the nearest 

integer. Consequently the row or column for 'total' is not necessarily the exact total. 



(b) the converse of that. I am not claiming that both halves 
of the model are correct, but only that they represent a 
reasonable starting point. Most stellar population synthesis 
calculations discount entirely multiplicity greater than two, 
and I believe the present model is somewhat better than 
that. 

The distribution over spectral type (second-last column 
of Table 5) is the same as in the last two rows ('multiple') in 
Table 1, apart from the fact that Table 1 contains the rms 
spread as well as the average. The only difference between 
the model in Table 5 and the 'OB down' model of Table 1 
is that multiples are now allowed. Probably the remaining 
discrepancies can be improved by varying the IMF a little. 
However there are at least two other ways of altering them: 
(i) 'convective overshooting', an effect modeled in my stellar- 
evolution code on the lines of Schroder et al. (1997) and Pols 
et al. (1997), can vary the length of time spent on the MS, 
differentially as a function of mass; and (ii) small changes, of 
order 50K, in the temperature boundaries of spectral types 
in Table 3 can send quite a lot of systems from say type A 
to type F. I prefer for the time being the rather considerable 
simplicity of the model used here, without trying to further 
improve the agreement. 

The selection effects implemented in Observe are very 
crude, but they seem to imply that in order to get approxi- 



mately the observed average multiplicity of 2.42 for O stars, 
we need to start with an average well in excess of 3.0. This 
average has to drop fairly rapidly to about 2 at later spectral 
types in order that the observed multiplicity should drop to 
~ 1.6, and even lower for M stars. The average multiplic- 
ity in the entire ensemble, before the Observe stage, was 
2.04 ± 0.03. After Observe, it dropped to 1.55 ± 0.02. In 
the HBMSC, and of course also the TOBMSC, triples and 
higher multiples constitute about 9% of systems, but in the 
RTBMSC they constitute 20%, which I suggest is a more 
realistic assessment of their frequency. 

In Table 6, the 'observed' periods of the HBMSC that 
are fairly long, in particular those over ~ 200 yr, are esti- 
mates that are based on Kepler's Law, a measured angular 
separation and parallax, and a mass based on spectral type. 
They are obviously very uncertain, but perhaps by not much 
worse than one bin (a factor of 10) either way. We see that 
the effect of selection, for OB stars, is to turn a largely uni- 
modal period distribution as seen in the RTBMSC section 
of Table 6 into a strongly bimodal distribution as seen in 
the TOBMSC and HBMSC sections. The 'missing' systems 
at periods of ~ 1 — 100 yr are presumably too wide (and 
broad-lined) for radial velocity curves, and too close for vi- 
sual resolution. 

In fact a slight bimodality can be seen even in the 
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RTBMSC period distribution for type O, which comes from 
the fact the the shorter periods are more probably from inner 
pairs via equation (18) than from outer pairs via equation 
(17). But the bimodality is much more marked when the 
selection effects are allowed for, and runs all the way from 
type O to type A, with a minimum at about 10 — 100 yrs as 
observed. 

The most obvious discrepancy between the HBMSC and 
TOBMSC sections of Table 6 is in the bottom left-hand 
corner, for types G/K/M and periods £s lyr. The fact that 
there are many in the TOBMSC and and few in the HBMSC 
is due to binary interaction. I have emphasised that this is 
not (yet) included in the Create code, and so the model 
is comfortable producing say a 0.01 yr binary containing a 
well-evolved red giant, that ought to have filled its Roche 
lobe some time ago. Such a system might now resemble a 
Leo (Gies et al. 2008) where a late B subgiant has a low- 
mass companion that is (probably) a white dwarf, in a 40.1 d 
circular orbit. This orbit was not known previously, and is 
not included in ET08. 204±27 such systems can be expected. 
Some will probably have merged into single stars, and others 
to some other part of the Table. 17 binaries or sub-binaries 
are known in the HBMSC to be semidetached, i.e. currently 
undergoing RLOF. A further handful, such as <j> Per, 7 Cas, 
9 Tuc and AY Cet, and now a Leo, are believed to be post- 
RLOF objects. But we can be confident that there are many 
more in this category that would be hard to recognise by any 
means. 

Leaving aside the bottom left-hand corner, we can see 
that there is reasonable agreement over the rest of the Table. 
I do not believe it would be helpful to try and quantify this 
by some x 2 statistic, or to overegg the pudding by trying 
to make the agreement even better with adjustments to, for 
instance, the period distributions (17) and (18). 

There appears to be another discrepancy: the total 
number of orbits in the HBMSC is 1979, and in the 
TOBMSC is 2315. But this is because there are several sys- 
tems or subsystems (in fact 384) in the HBMSC that do not 
have a measured, or even estimated, period. These are for 
example the 'astrometric accelerators' of Makarov & Kaplan 
(2005), radial velocity variables where probably the motion 
is orbital but no orbit has been determined (Nordstrom & 
Andersen 1985) and about 40 Ba stars for which a white- 
dwarf companion is expected and yet no orbit has yet been 
measured. 

One of the more important uncertainties in the model is 
the mass-ratio distribution (19). It is unlikely that a single 
value of the exponent a will apply independently of m, P 
and hierarchical depth I. It would be desirable to check this 
much more rigorously, by using observed mass ratios. How- 
ever in observational analysis the mass ratio is one of the 
more difficult quantities to measure. A spectroscopic binary 
has to be double-lined, with the secondary preferably not 
very faint so that its amplitude can be well determined. 
This means that mostly only systems with Q fairly close 
to unity are measurable, and so the regime of Q ^ 0.5 is not 
well tested. The TOBMSC model of Tables 5 and 6 produced 
198 ±21 double-lined and 854 ±56 single-lined spectroscopic 
systems or subsystems. In the HBMSC 638 systems or sub- 
systems are SB1 or SB2. The discrepancy is fairly large, 
but is quite probably because all (theoretical) systems or 
subsystems with P ^ 50 yrs and large enough radial veloc- 



ity amplitude (Table 3) are recognised as SB, whereas it 
is likely that many real systems in the period range of say 
10 — 50 yrs have not been recognised. I am not claiming that 
a global choice of a — 0.8 is correct, but only that it works 
reasonably well to produce an ensemble like that observed. 

The TOBMSC contained 85 ± 12 eclipsing binaries, 
against 130 in the HBMSC. The RTBMSC (Table 6) con- 
tained as many as 32 very short-period (P < 0.001 yr) bina- 
ries, most of which were eclipsing; but these were all double- 
M-dwarf sitfe-binaries of systems whose primary had the type 
listed, and all but two on average were undetectable in the 
Theoretical Observatory. The HBMSC contains two systems 
at such short periods, both being subsystems in fact. The 
well-known M-dwarf eclipsing sub-binary in the sextuple 
q Gem is almost in the same category, but its slightly longer 
period puts it in the 0.001 — 0.01 yr bin. 

Barium-rich (Ba) stars are thought (McClure 1983, Bof- 
fin &Jorissen 1988) to be binaries where the primary is now 
a white dwarf, and where the secondary, a G/K giant, has 
been contaminated from s-processed material in the enve- 
lope of the white dwarf's precursor when it was near the 
tip of the Asymptotic Giant Branch (AGB). There are 53 
Ba stars in the HBMSC: 11 of them have known SB orbits, 
with periods ranging from 285 to 6500 days; and 3 of them 
have known white dwarf companions, notwithstanding the 
difficulty of recognising a white dwarf companion whose lu- 
minosity will usually be one part in 10 3 or 10 4 of the red 
giant. In the TOBMSC, all systems containg a white dwarf, 
a G/K giant, and with period in the range 100 — 10000 d 
were identified as potential Ba stars. There were 37±15, as 
compared with the 53 in the HBMSC. But many Ba stars 
in the HBMSC are only mildly enriched in Ba, and if this 
can be achieved in wider binaries with periods up to 20000 d 
then there is no discrepancy (56 ± 25). 

I emphasise again that the HBMSC only includes sys- 
tems more massive than about 1 Mq in total. Coincidentally, 
but conveniently, this is also the set of stars that are capa- 
ble of evolving significantly in the course of the Galaxy's 
lifetime. Thus this model may be quite adequate for explor- 
ing the evolution of an ensemble of multiples (including by 
definition singles) on the scale of a galaxy, or a galactic clus- 
ter. In the context of a cluster, we can be reasonably sure 
that the wider systems will be broken up quite quickly, pro- 
vided that the cluster remains rich for a reasonably long 
time. However, current thinking is that all stars form in rich 
or fairly rich SFRs. This must therefore include the rather 
highly- multiple, yet also quite wide, O star multiples of the 
HBMSC, which may indeed be the result of (a) gravitational 
diffusion inwards of the most massive stars of a cluster, and 
(b) dissipation of the cluster because of the ejection of signif- 
icant mass by, for instance, supernova explosions. I believe 
that one of the more rigorous tests of star-formation models 
will be that they produce the right spectrum of multiplicity, 
as well as the right spectrum of masses and orbital periods. 

A striking result of rather recent analyses (Tokovinin 
et al. 2006, Pribulla &Rucinski 2006), confirmed by ET08, 
is that short-period binaries (P < 3d) are very commonly 
in triples, and arguably are exclusively in triples (or higher 
multiples); and this then implies that tripleness is neces- 
sary for the production of a short-period binary. The nec- 
essary mechanism may be the combination of Kozai cycles 
and tidal friction: KCTF (Mazeh & Shaham 1979, Kiseleva 
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ct al. 1998, Eggleton & Kisseleva-Eggleton 2006, Fabrycky 
& Tremaine 2007). It may be that the process of star for- 
mation does not directly produce binaries with periods of 
less than say 0.1 — 1 yr, but it produces enough triples, with 
suitably oblique orbits, that these short periods can be gen- 
erated by the KCTF mechanism, in the course of 10 3 — 10 9 yr 
depending on the initial parameters. Once an orbital period 
has been reduced to ~ 2 — 3d by KCTF, those systems 
which contain late- type dwarfs (F/G/K/M) may evolve to 
still shorter periods by the quite different process of mag- 
netic braking and (again) tidal friction, or MBTF. This can 
reduce the period to the point where a contact binary is 
formed, on a timescale that may be roughly 10 7 — 10 10 yr. 

The reader will note that there is something of a con- 
tradiction between 

(a) offering a formulation that is intended to give roughly 
the right distribution of periods at age zero, 

(b) generating arbitrarily short periods via equations (17) 
and (18), although with the distribution truncated when the 
period is so short that the stars would be touching on the 
ZAMS, and 

(c) hypothesising that the shortest periods ( Ss 3 d) are in 
fact due to the KCTF and MBTF mechanisms, and not to 
formation at such short periods. 

I would hope, in a future analysis, to investigate whether 
by considering the interactive evolution of triples hypothe- 
sis (c) can be supported, but until then it seems desirable 
to produce about the right frequency of such short-period 
systems ab initio. 

Part of the above apparent contradiction is because 'age 
zero' is rather hard to define. Returning to the brief discus- 
sion of star formation in the introduction, the bright stars in 
the Solar neighbourhood presumably formed in SFRs, and 
then escaped as the SFRs evaporated; but they may have 
been subject to dynamical interaction within their parent 
SFRs for at least several million years, so that their multi- 
plicities and periods may have been modified in that inter- 
val. 'Age zero' can be surprisingly well defined for individual 
stars when their nuclear evolution alone is being discussed. 
But 'age zero' is much harder to define when binaries and 
higher multiples are being considered. 

The purpose of this paper has been to produce a 
model for the multiplicity statistics of an observed com- 
plete magnitude-limited sample of substantial size. Such a 
model must include selection effects, but even the best such 
model will not be definitive. In particular, one can add to al- 
most every observed system a T dwarf that would hardly be 
recognisable by current technology, if located at a suitable 
distance from the primary component. We can only attempt 
to find something like a lower limit to multiplicity; but I be- 
lieve that the model proposed here is something like a lower 
limit. 
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